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The cutoff method, which cuts off the values of a function less than a given number, 
is studied for the numerical computation of nonnegative solutions of parabolic partial dif- 
ferential equations. A convergence analysis is given for a broad class of finite difference 
methods combined with cutoff for linear parabolic equations. Two applications are in- 
vestigated, linear anisotropic diffusion problems satisfying the setting of the convergence 
analysis and nonlinear lubrication-type equations for which it is unclear if the convergence 
analysis applies. The numerical results are shown to be consistent with the theory and in 
good agreement with existing results in the literature. The convergence analysis and appli- 
cations demonstrate that the cutoff method is an effective tool for use in the computation 
of nonnegative solutions. Cutoff can also be used with other discretization methods such 
as collocation, finite volume, finite element, and spectral methods and for the computation 
of positive solutions. 
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1 Introduction 

Many physical phenomena involve variables such as the density and concentration of a material that 
take only nonnegative vahics. The mathematical reflection of this property is that the partial dif- 
ferential equations (PDEs) modeling the phenomena admit nonnegative solutions. For the numerical 
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solution of those PDEs, it is crucial that numerical schemes preserve the solution nonnegativity and 
produce physically meaningful numerical solutions. 

A closely related concept is the maximum principle. Preserving the maximum principle is equivalent 
to preserving the solution nonnegativity for linear problems [50] but generally the former is more 
difficult than the latter. It is known (e.g., see that a conventional numerical method, such as a 
finite difference (FD), finite volume (FV), or finite element (FE) method, generally does not preserve 
the maximum principle and can produce negative undershoot in the solution for diffusion problems, 
especially those with heterogeneous anisotropic diffusion coefficients. Considerable effort has been 
made on developing numerical schemes satisfying the maximum principle; e.g., see [U [71 \T0\ \TT\ [28l 
[29l[3ll[3lll6l[171[5ll[52| for steady-state isotropic diffusion problems and[l5l[22l[23l[271[30l[36l[38l[39l 
Uni HU [321 US] for steady-state anisotropic diffusion problems. Particularly, Ciarlet and Raviart [TT] 
show that the linear FE approximation of an elliptic isotropic diffusion problem satisfies the maximum 
principle when the mesh is simplicial and nonobtuse. The result is generalized to anisotropic diffusion 
problems by Huang and his coworkers in |27 | [36 l E]. On the other hand, less progress has been made 
for time dependent problems. For example, Fujii [21] shows that the linear FE approximation of the 
heat equation satisfies the maximum principle when the mesh is simplicial and acute and the time step 
size is bounded above by a bound proportional to the squared maximal element size and below by a 
bound proportional to the squared minimal element in-diameter. He also shows that when the lumped 
mass matrix is used, the maximum principle holds without requiring the time step size to be bounded 
from below. Fujii's results are extended to more general isotropic diffusion problems by Farago et 
al. [T71 |T9l |20] and to anisotropic diffusion problems by Li and Huang [37]. Farago and Horvath 
|18j study the relations between the maximum principle, nonnegativity preservation, and maximum 
norm contractivity for linear parabolic equations and their finite difference and Galerkin finite element 
discretizations. Thomee and Wahlbin [50j consider more general parabolic PDEs and show that the 
maximum principle cannot hold for the conventional semidiscrete FE problem and Fujii's conditions 
on the mesh are essentially sharp for the lumped mass matrix. Le Potier |32[ [55] (also see Lipnikov et 
al. [39] ) proposes two nonlinear FV schemes for linear anisotropic diffusion problems and shows that 
they are second order in space and satisfy monotonicity [32j or the maximum principle [33j. 

We consider initial-boundary value problems of parabolic PDEs which are well posed and have a 
unique, nonnegative solution. Instead of employing a maximum-principle preserving scheme, we pro- 
pose to use the cutoff method that cuts off the negative values in the computed solution at each time 
step and then continues the time integration with the corrected solution; see Fig. [T] for a sketch of the 
solution procedure. The cutoff method shares a similar idea with many projection and correction meth- 
ods. For example, projection methods [HlllHj are an effective means of enforcing the incompressibility 
condition in the numerical solution of time dependent incompressible fluid flow problems. Projection 
is also used by some researchers in the numerical solution of Hamiltonian systems to preserve the 
energy; e.g., see [26]. A EE/implicit Euler discretization of a Cahn-Hilliard equation with logarithmic 
nonlinearity is considered in [13] and iterative methods to solve the resulting nonlinear systems of 
equations are developed to keep the solution within its physically relevant range. Solution compres- 
sion is used in the design of maximum-principle-satisfying high order schemes for scalar conservation 
laws by Zhang and Shu [53] . It should be emphasized that the cutoff method provides several advan- 
tages over many positivity-preserving or maximum-principle-preserving schemes. Its implementation 
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is simple and requires no significant changes in the existing code. Moreover, the cutoff procedure does 
not impose any constraint on the mesh and time step which maximum-principle-preserving schemes 
often impose. However, unlike projection methods, there do not seem to exist published theoretical 
or numerical studies on the cutoff method. 

The objective of this paper is to provide a theoretical and numerical study of the cutoff method 
for use in the computation of nonnegative solutions of parabolic PDEs. We shall first prove that a 
broad class of finite difference methods are convergent when they are incorporated with the cutoff 
method for linear parabolic PDEs. Two applications are then investigated. The first is a linear 
anisotropic diffusion problems which satisfies the setting of the convergence analysis. It is known 
that a conventional finite difference or finite element method with a uniform or quasi-uniform mesh 
does not preserve the maximum principle (nor the solution nonnegativity) and typically produces 
negative undershoot in the computed solution. The cutoff method removes the unphysical negative 
values in the solution while keeping the same convergence order of the underlying scheme. The second 
application is nonlinear lubrication- type equations. It is unclear if the convergence analysis (which 



can be extended to some nonlinear equations; see Remark 2.1) applies to those equations. A distinct 
feature of lubrication-type equations is that a positive solution can develop a finite time singularity 
of the form n — t- and become identically zero on some spatial interval for a period of time. A 
conventional finite difference or finite element method can produce negative values in the computed 
solution during this development of singularity. Negative solution values are not only unphysical 
but also cease the computation for typical situations where the nonlinear diffusion coefficients are 
undefined for negative solution values. The application of the cutoff method avoids this difficulty. 
Moreover, it is shown that no regularization is needed when cutoff is used. Numerical results are 
shown to be in good agreement with existing ones in the literature. 

The outline of the paper is as follows. The cutoff method is described and some of its properties 
are given in Section [2j The convergence of a broad class of finite difference methods is proved also 
in Section [2] when they are incorporated with the cutoff method for the computation of nonnegative 
solutions of linear parabolic PDEs. Application of finite difference methods with cutoff to a linear 
anisotropic diffusion problem and a lubrication- type PDE is studied in Sections [3] and [4j respectively. 
Finally, Section [5] contains conclusions and further comments. 

2 The cutoff method and convergence analysis 

In this section, we first describe the cutoff method and some of its properties. We then give a 
convergence analysis of the method for a class of finite difference methods applied to linear parabolic 
problems. Possible generalizations are also discussed. 

2.1 The cutoff method 

Consider a continuous function / = f{x) defined on a domain C (d > 1). For any given cutoff 
parameter 5 € M, we define the (5-cutoff function as 

_i_ , ^ I / for f(x) > 6 

fli^) = { : ' , ' - ea (1) 

for j{x) < 
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Figure 1: An illustration of the time integration with cutoff (the dashed line). Here, C/" is the numer- 
ical solution, (C/")^ is the corrected numerical solution and u{t) is the exact solution of the 
IB VP starting from at t = to- 



Note that f^{x) = fQ{x) is the nonnegative part of /, i.e., 

2 [O, for f{x) < 

The following three lemmas give some properties of cutoff functions. Although these properties are 
easy to prove (and almost obvious), they play a central role in the convergence analysis for a class of 
finite difference methods in the next subsection. 

Lemma 2.1. For a given function f and any nonnegative continuous function u defined on 0,, we 
have 

\f+{x)-u{x)\<\f{x)-u{x)\, yxen. (3) 

Proof. Obviously, ([s]) holds when f{x) > 0. When f{x) < 0, we have f^{x) = 0. By assumption, 
u{x) > 0. Thus, we have 



\f^{x) - u{x)\ = \u{x)\ < \f{x)\ + \u{x)\ = \fix) - u{x)\. 



n 



Lemma 2.2. For a given function f and any nonnegative continuous function u defined on 0,, we 
have 

\f+{x)-f{x)\<\u{x)-f{x)\, yxen. (4) 

Proof. The proof is similar to that of Lemma |2.1[ □ 

Lemma 2.3. Given a cutoff parameter 6 >0, for any function f and any nonnegative continuous 
function u defined on Cl, we have 



\f+{x)-f+ix)\<6, yxen 

\f+{x) - u{x)\ < \f{x) - u{x)\ + 6, yxen. 



(5) 
(6) 
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Proof. Inequality ([5j) follows from the definitions of and f^. Inequality ([6j) follows from the 
triangle inequality \fg~{x) — u{x)\ < \f~^{x) — u{x) \ + |//(£c) — f^{x)\ and Lemma 



2.1 



2.2 Error analysis for a linear IBVP problem 

We now use the results in the previous subsection to analyze the convergence of a class of finite 
difference schemes for a general initial-boundary value problem (IBVP) in the form 



^ = L{u), in J^x(io,T] 
u = g, on 

u = u^, on Q, X {t = to} 



(7) 



where C M*^ (d > 1) is the physical domain, L is a linear elliptic differential operator and g and 
are given sufficiently smooth functions. We assume that the IBVP is well posed and admits a unique 
nonnegative continuous solution, u = u{t, x). We also assume that L and g do not contain t explicitly 
for notational simplicity. 

We consider a class of finite difference schemes for ([T]) in the matrix form as 



(8) 



where Bq and Bi are matrices independent of n and C/" is an approximation of the exact solution at 
t = tn, i.e., = {Uj'}, = {Uj}, Uj- ^ Uj = u{tn,Xj) for each mesh vertex Xj. We assume that 
([s]) satisfies 



1 II — -''■ii 
Bi^BoW < (l + KAt), 



1 

At 



lU 



71+1 



Bqu'' - F'- 



du 
di 



L{u), as At{h)^0 



(9) 
(10) 

(11) 



where At and h are the maximal time step size and the maximal element size, respectively, K and 
Ki are positive constants, and || • || is a proper matrix norm. Condition ^ requires that ([s]) produce 



a bounded solution while (10) and (11) are the stability and consistency conditions, respectively. The 



local truncation error of this scheme is defined as 



(12) 



Assume that scheme ([s]) is {p, q)-order for some positive integers p and q. Then, there exists a constant 
Cite{u) (depending only on the exact solution) such that 



|t"|| <Qte(n)At(AtP + /i^). 



(13) 



It is remarked that there exist many schemes satisfying assumptions ^ - (11). Examples include 
those employing central finite differences for spatial discretization and the ^-method for temporal 
discretization; e.g., see Morton and Mayers 
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Theorem 2.1. Assume that IBVP ^ is well posed and admits a unique nonnegative continuous 
exact solution u = u{t,x). We also assume that scheme satisfies ^ - (11) and (13). Then, the 
■ for the cutoff solution procedure shown in Fig. [7] with scheme is hounded by 



error 



||(C/")+ - ^/"ll < ||[/" - < ^ e^*" Cite{u) {AtP + h"). 

K 



(14) 



Proof. Let e" = C/" — ti". Notice that the cutoff solution procedure shown in Fig. [T] with ([s]) 
satisfies 

= Bo(t/")+ + i^". (15) 



Combining this with (12), we obtain 

i?ie"+i = i?o((t/") 

It follows that 



Combining this with \Sm and ( 10 ) leads to 



Lemma 2.1 implies 



Thus, we get 



e"+i < {l + KM) II ([/")+ -u" II ||r" 

11(^7")+ -u'^ II < ||?7"-u"||. 
||e"+i|| < (1 + KM) ||e"|| + Ki ||r"||. 



Then it is standard to show that (14) follows from (^13| and (18). 



(16) 

(17) 

(18) 

□ 



Remark 2.1. From the above proof we can see that the key condition is the convergence of the 



original scheme (without cutoff). If it (without cutoff) is convergent, using Lemma 2.1 we can readily 
show that the scheme with cutoff is also convergent. This observation implies that the convergence 
analysis of the cutoff method can be extended to more general linear or nonlinear IB VPs. □ 



When a strictly positive solution is desired, we can replace ([/")+ with (t/")t for some positive 



constant 6. The following theorem can be proven using Lemma [2. 3| in a similar way as for Theorem 2.1 



Theorem 2.2. Assume that IBVP ^ is well posed and admits a unique nonnegative continuous 
exact solution u = u{t,x). We also assume that scheme satisfies ^ - (11) and and (IS). For the 
cutoff solution procedure shown in Fig. [7] with scheme ^) and with ([/")+ being replaced with {U"-)^ 
for some positive 6, the error is bounded by 



u 



< — e 
- K 



Cueiu) [At^ + /i^) + ( e^*" + 1 U. 



(19) 



Remark 2.2. Inequality (19) shows that 6 should be chosen proportional to At for the error 



bound to stay bounded as At — )• 0. Ideally, 5 should be at the same level as the local truncation error, 
i.e., 

5 = 0(At(AtP + /!«)). (20) 



This way, the terms on the right-hand side of (19) have the same convergence order as At(/i) — t- 0. n 
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3 An anisotropic diffusion problem 



In this section we present numerical results obtained for a 2D linear anisotropic diffusion problem by 
the cutoff method described in the previous section. The problem takes the form of IB VP Q with 

L{u) = V-iBVu) + f, n=[0,l]x[0,l], ©=(^^43^ /qq^ ^ , (21) 

and /, u", and g are chosen such that the exact solution of the IB VP is given by 

u = ^exp(-t)(tanh(-15(x - y)) + 1). (22) 

The problem satisfies the maximum principle and the solution stays between and 1. 

It is worth pointing out that unlike lubrication-type equations which we shall consider in the next 
section, for this problem the computation can continue when negative values occur in the computed 
solution. From this point of view, it is unnecessary to remove negative values at each time step. 
Nevertheless, this problem satisfies the setting of the convergence analysis in the previous section and 
provides a good example for verifying the theory and testing the effectiveness and accuracy of the 
cutoff method. 

We use Cartesian grids of size J x J for the physical domain, central finite differences for spatial 
discretization, and a third-order singly diagonally implicit Runge-Kutta (SDIRK) method [U [H] for 
temporal discretization of the underlying PDE. The discretization is standard and can be shown to 
be convergent with order (3, 2) (third order in time and second order in space). It is also known (and 
is shown below) that this standard scheme does not preserve the maximum principle and produces 
spurious undershoot and overshoot in the numerical solution. The numerical results presented below 
are obtained with a fixed time step size At = 10~^, which was found to be small enough so that the 
temporal discretization error is ignorable compared to the spatial discretization error. 

Fig. [2] shows the surface plot of a numerical solution (after cutoff) at t = 1. The contours of the 
numerical solutions before and after cutoff are shown in Fig. [s] (a) and (b). One may notice that the 
underlying finite difference scheme does not preserve the nonnegativity and the numerical solution at 
t = 1 contains negative values (with C/^j„ = —1.073 x 10"'^ where is the last time step) before 
cutoff. 

The error, IK?/^)''' — u^\\l^(q), and the maximal undershoot, —Umin, are shown in Fig. jZ] as 
functions of the number of subintervals in x (or y) direction. We can see that the error decreases 
at a rate of 0{J~'^) (second order in terms of element size) as J increases. This is consistent with the 
theoretical prediction given in Theorem |2.1[ On the other hand, the maximal undershoot, which is 
equal to the cutoff error, i.e., —Umin = — (t^^)'*'||oo) decreases at a faster rate. This is somewhat 
surprising since we expect the cutoff error to be at the level of the local truncation error, which is 
second order in space. 

Recall that the cutoff strategy can be used to preserve the positivity of the solution. For example. 



when the cutoff parameter is taken 6 = At(Ax)^, we have {U^)~^ > At(Aa;)^. Theorem 2.2 implies 
that the error, ||(f7''^)^ — decreases at a rate of second order in space. For the current 

example, the numerical results obtained with this cutoff parameter are indistinguishable from those 
shown in Figs. [3]and[4j For this reason and to save space, we omit those results here. 



7 



It is worth mentioning that computations were also performed for a similar problem with a convec- 
tion term, 

L{u) = V • (D Vm) - 5 • Vu + /, (23) 
where D is given in py, b = [1000,1000]^, and the functions /, g, and are chosen such that the 



exact solution is given by (22). The obtained results (not shown) are very similar to those shown in 
Figs. [3]and[4l 








Figure 2: Example (21 ). The surface plot of a numerical solution (after cutoff) at t = 1 obtained with 
an 80 X 80 Cartesian grid. 



4 Application to lubrication-type equations 

In this section we consider the application of the cutoff method to a lubrication-type equation, which 
was first derived from lubrication theory by Greenspan |24j for describing the movement of thin 
viscous films and spreading droplets. Lubrication theory consists of a depth-averaged equation of 
mass conservation and a simplified form of the Navier-Stokes equations that is appropriate for a thin 
layer of very viscous fiuid. We consider a lubrication-type equation in the general form 



ut + V ■{fiu)VAu) = 0, f{u) 



u 



n G [0, cxd) 



(24) 



where u is the thickness of the viscous droplet and n is a physical parameter which has different values 
for different boundary conditions on the liquid solid interface. Lubrication- type equations also appear 
in several other applications, including a thin neck model in the Hele-Shaw cell |12j, Cahn-Hilliard 
models with degenerate mobility jT6], and problems in population dynamics [35] and plasticity [25] . 
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(a): U- 



(b): {U^) + 




It is noted that the solution u must stay nonnegative to be physicaUy meaningful in all of those 
applications. 



It is known theoretically [21 E] that in one dimension, (24) preserves the positivity of the solution 
for n > 3.5 and has a nonnegative weak solution in a sense of distributions for 3/8 < n < 3 (and in a 
weaker sense for < n < 3/8). Such a weak solution can be obtained [2] as the limit as e — t- of the 
solution of the regularized problem 



e^4 



e\4 ■ 



(25) 



where e > is the regularization parameter. Moreover, numerical computation (such as see [3ll3l[5]) 
shows that for small values of n > 0, a positive solution can first develop a finite time singularity of 
the form u — )• at some point (which physically describes the rupture of the liquid film) , then becomes 
identically zero on an interval of time dependent length for a period of time, becomes positive again 
at a later time, and eventually decays to the mean of the initial solution. It is challenging to simulate 
this singularity development since conventional numerical methods do not preserve the nonnegativity 
of the solution and requires a huge number of mesh points to provide the necessary resolution to 
avoid spurious, negative undershoot. Moreover, when n is an even denominator rational number, 
f{u) is undefined for negative u. In this situation, the computation typically stops around the initial 
formation of the singularity when negative values start to appear in the numerical solution. For this 



reason, the simulation of the singularity development in (24) is a good test for the cutoff method 



although it is unclear if the convergence analysis described in § 2.2 applies to (24). 



Great effort has been made in the numerical solution of (24). In addition to the above mentioned 
references [3 HI E]; Zhornitskaya and Bertozzi [S3] propose a positivity preserving finite difference 



scheme for the regularized equation (25). Russell et al. [31] solve the same regularized equation using 



a moving collocation method. Sun et al. [48J solve (24) in two dimensions using an adaptive finite 
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Figure 4: Example (21). The L? error, \\{U^) 



L2(^Q^, and the maximal undershoot, —Umin, are 



shown as functions of the number of subintervals in x (or y) direction. 



element method and show that proper mesh adaptation can provide accurate resolution and there is 



no need to use regularization in the numerical simulation of the singularity development in ( 24 ) . 
We first consider an IBVP of the one-dimensional lubrication-type equation as 



Ut + {f{u)Uxxx)x = 0, X G r2 EE (-1, 1), f{u) 
U{t, ±1) = 0, Uxxx{t, ±1) = 0, 

u{0, x) = 0.8 - cos(7rx) + 0.25 cos(27rx). 



1 

U2 . 



(26) 



Notice that the initial condition is strictly positive (with minimum value 0.05 and mean value 0.8) and 
f{u) is undefined for negative u. We use a uniform mesh of J cells for il, central finite differences for 
spatial discretization, a third-order SDIRK method for the temporal discretization of the underlying 
PDE. (A small fixed time step At = 10~^ is used in our computation.) This discretization is basically 
the same as that for the example in the previous section except that a special treatment is needed 
for the diffusion coefficient for the current problem since it is nonlinear. Recall that the scheme 
does not preserve the nonnegativity of the solution. Thus, iterates can have negative values at some 
point during the Newton iterative solution of the nonlinear algebraic equations resulting from the 
SDIRK/finite difference discretization of the underlying PDE. Once this occurs, the computation 
stops because f{u) is undefined for those values. One way to avoid this difficulty is to use f{{U"^^)~^) 
instead of f{U"'~^^) in the scheme. However, the nonsmoothness of {U"'~^^)~^ can cause difficulty in the 
convergence of the Newton iteration. We use here the "lagged diffusivity" method, i.e., the diffusion 
coefficient is computed using the (corrected) numerical solution at the previous time step. An iterative 
method based on the lagged diffusivity is used and proved to be convergent by Dobson and Vogel [H] 
for a total variation denoising problem which also has a nonlinear diffusion coefficient. 

Fig. [sjshows a numerical solution at various time instants obtained with the cutoff method (without 
using regularization for f{u)) on a relatively coarse uniform mesh of 129 points. The result is almost 
identical (in the eyeball norm) to those in [3l SH]. To further verify the result, the computation is 
done with a uniform mesh of 1001 points and the obtained solution is indistinguishable from the one 
shown in Fig. [5] in the plotting precision. 
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A uniform mesh of 1001 points is used for the numerical study of the singularity development of 
the solution. Fig. [6] shows the close views of the solution during the development. It can be seen that 
the singularity is developed at around t = 7.30 x 10"^, which is consistent with existing numerical 
simulations, including those in [53] with a positivity-preserving scheme. Then, the solution becomes 
identically zero on an interval of time dependent length for a period of time. The length of the interval 
increases from zero, attains its maximum (about 0.12), and decreases to zero. Afterward, the solution 
becomes positive again, and eventually decays to the mean of the initial solution. It is known (e.g., 
see [3]) that the development of the singularity is characterized by the third order derivative of the 
solution becoming infinite at certain points. This can be seen in Fig. [7]where the third order derivative 
of the numerical solution is shown at various time instants. In particular, the third order derivative 
becomes discontinuous when the singularity occurs, then has jump discontinuities, and then becomes 
continuous again when the solution is positive. 



As we have seen in the above, there is no need to use any regularization for f{u) (cf. (25)) with the 
cutoff method. On the other hand, it is interesting to see the effects of the regularization of f{u) on 
the solution. We first point out that the regularization does not guarantee the nonnegativity of the 
solution for the central finite difference discretization on a uniform mesh. Thus, we also need to use 
the cutoff method for the regularized equation (25). A solution obtained for (25) (with e = 10~^^) 



with the cutoff method is shown in Fig. ^ By a direct comparison with Fig. [Tj one can see that there 
is only a slight difference (at the level of 0(10"^)) in the onset time and disappearance time of the 
singularity. Moreover, the solution of the regularized problem does not become identically zero on 
an interval for a period of time. Instead, there is a bump at the central part of the interval and the 
height of the bump maintains constant almost for the whole appearance period of singularity. Once 
again, the length of the interval first increases from zero, reaches its maximum, and decreases to zero 
when the solution becomes positive again. 

Fig. [9] shows comparison of numerical solutions at onset of singularity and at t = 0.001. It can be 
seen that the regularization changes the onset pattern: the solution to the non-regularized equation 
touches the x-axis at one point whereas those to the regularized one touch at two points simultaneously. 
Moreover, the length of the touching interval and the height of the bump at onset depend on the 
regularization parameter e. The smaller e, the smaller the touching interval and the lower the bump. 
This suggests — t- n as e — t- 0. It is interesting to point out that all solutions, for regularization or 
non-regularization, have almost the same length of the touching interval at t = 0.001. 

Finally, we present some numerical results for an IBVP of the two-dimensional lubrication-type 
equation as 



1 



ut + V ■{f{u)VAu) = 0, {x,y)en = [-l,l]x[-l,l], /(n)=n2, 
u = 0, ^ = 0, {x,y)edn 

u{0, X, y) = (0.8 - cos(7rx) + 0.25 cos(27rx))(0.8 - cos(7ry) + 0.25 cos(27ry)), (x, y) G 



Fig. 10 (a) shows a numerical solution at t = 0.001 obtained with a uniform mesh of size 81 x 81 while 



its profiles along the diagonal line x = y at various times are shown in Fig. 10 (b). One can see that 
the development of singularity has a similar behavior as in one dimension (cf. Fig. [s]). One may also 
see that the onset time for the current example is slightly earlier. This is because the current initial 
solution has a smaller minimum value 0.025. The contours of the numerical solutions at i = 0.001 
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before and after cutoff are shown in Fig. 11, The undershoot (with Umin — ~ 

3.23 X 10""^) is visible 

in the solution before cutoff. These results demonstrate that the cutoff method can be used for the 
simulation of singularity development in the two-dimensional lubrication-type problem without need 
of any type of regularization. 




0.5 1 



Figure 5: The one-dimensional lubrication-type problem. The solution obtained with 129 uniform grid 
points is shown at various time instants. 



5 Conclusions and further comments 

In the previous section we have studied the cutoff method for the numerical computation of nonneg- 
ative solutions of parabolic partial differential equations. Several properties of the cutoff method are 



given in Lemmas 2.1 - 2.3 Convergence of a class of finite difference methods is proved (Theorems 2.1 
and 2.2 ) when they are incorporated with the cutoff method for linear parabolic PDEs. The method is 
investigated for two applications, linear anisotropic diffusion problems and nonlinear lubrication-type 
PDEs. The numerical results are consistent with theoretical predictions and in good agreement with 
existing results in the literature. 

We have considered finite difference methods in this work. But it is worth pointing out that the 
cutoff method can also be used with other discretization methods, e.g. collocation, finite volume. 



finite element, or spectral methods. As an example, we show in Figs. 12 and 13 results obtained with 
the standard linear finite element method on Delaunay meshes for anisotropic diffusion problem (21). 
They are comparable with those in Figs.[3]and|4j Theoretical analysis for finite element methods with 
the cutoff strategy is currently underway. 
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Figure 6: The one-dimensional lubrication-type problem. The close views of the numerical solution are 
shown at various time instants during the singularity development. The numerical solution 
is obtained by the cutoff method (without using regularization for f{u)) on a uniform mesh 
of 1001 points with At = 10"^. 
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